%Read MR Solutions dynamic MRS data
clear all;  close all;

global sp

sw=3330;    ivt=24;   ppmL=160;   ppmH=195;   ave=1;
filt=1;     w=0.9; % filter weight-not used
ft_pts=512;
zerofill=0; % number of zeros to add at end
rank_dim=5 % degree of noise reduction. Set to the number of components, usually 2-5

% Select MRD Data Files
[File,Path,Filterindex] = uigetfile('*.MRD','Select MRD File','Multiselect','on'); 
cd(Path)
[rawdata,dim,par] = Get_mrd_3D1(File,'seq','seq');

%% FFT dim 2 to spectrum
pts=par.NO_SAMPLES;
ary=par.EXPERIMENT_ARRAY;
offset=par.freq_offset;

%% Tukeywin filtering -disabled as Signal Control box is missing
%if filt==1;
 %   f=tukeywin(pts*2,w);
  %  for i=1:ary
   %     F(i,:)=squeeze(f(pts+1:end));
    %end
    %data=data.*F;
%end


tmp=zeros(ary,zerofill);
data = horzcat(rawdata,tmp); %add zeros at end
data=data(:,1:512);
ft_pts=512;
sp=fft(data,[],2);


max_val = max(sp(:))
[max_x max_y]=ind2sub(size(sp),find(sp==max_val))





entropy_min=ACMEentropy_fun([0,0],sp(max_x,:),max_x);
for psi=1:360
    entropy_value=ACMEentropy_fun([psi,0],sp(max_x,:),0);
    if entropy_value<entropy_min;
        entropy_min=entropy_value;
        theta=psi;
    end    
    
end 

max_ppm=(offset-11498+sw/2)/32.1468+165.0;
min_ppm=(offset-11498-sw/2)/32.1468+165.0;

ppm = linspace(max_ppm,min_ppm,(ft_pts+zerofill));  % for 3T


for i=1:ary

    sp2=fftshift(squeeze(sp(i,:)));
    sp2=real(fft(ifft(sp2)* exp((-1)^(0.5)*theta*pi/180)));
    sp3(:,i)=sp2;

end
sp=flipud(sp3);
[U,S,V]=svds(sp,rank_dim);
sp=U*S*V';

%%background subtraction


sp_base=sp(mean(round(.66*ft_pts):round(.7*ft_pts)),: );
sp_base=repmat(sp_base,ft_pts,1);
sp=sp-sp_base;

% %Display dynamic spectra



figure;     hold on;    hidden on;
L=max(find(ppm>ppmL));  H=min(find(ppm<ppmH));
for i=1:(ary/2-ivt)
    j=H:L;
 
    x(H:L)=i*2;
    plot3(x(j),ppm(j),sp(j,i*2+ivt),'color','b')
end
hold off

smax=max(max(sp));
smin=min(min(sp));
xlim([0 (ary-ivt)]); ylim([ppmL ppmH]); zlim([smin smax]);
view([-150 50]);
xlabel('Time (s)'); ylabel('Chemical shift (ppm)'); zlabel('Signal intensity');

% %% Averaging
% for i=1:((ary-ivt)/ave)
%     sp_ave(:,i)=sum(squeeze(sp(:,(ivt+(i-1)*ave+1):(ivt+i*ave))),2)./ave;
% end
% 
% figure;     hold on;    hidden on;
% L=max(find(ppm>ppmL));  H=min(find(ppm<ppmH));
% for i=1:((ary-ivt)/ave)
%     j=H:L;
%  
%     x(H:L)=i*ave;
%     plot3(x(j),ppm(j),sp_ave(j,i),'color','b')
% end
% hold off
% smax=max(max(sp));
% xlim([0 (ary-ivt)]); ylim([ppmL ppmH]); zlim([smin smax]);
% view([-65 65]);
% xlabel('Time (s)'); ylabel('Chemical shift (ppm)'); zlabel('Signal intensity');
% 
% figure;     hold on;   
% xlim([0 (ary-ivt)]); ylim([ppmL ppmH]); zlim([smin smax]);
% xlabel('Time (s)'); ylabel('Chemical shift (ppm)');
% contour(sp,50);
% hold off;

% write to file
% file_name=strcat(Path,File(1,1:6),'spc');
% dlmwrite(file_name,sp,'delimiter','\t','precision',3)

% %% colormap
% 
figure;     hold on;    hidden on;
 c=linspace(10,1500,100);
contour(ppm(round(.3*ft_pts):round(.68*ft_pts)),linspace(1,240,240),sp(round(.3*ft_pts):round(.68*ft_pts),:)',c)
set(gca,'XDir','Reverse')


%%plotting metabolites
[M,I] = max(sp(220:230,50));
I_mal=I+219;
[M,I] = max(sp(250:260,50));
I_fum=I+249;



% [M,I] = max(sp(320:335,60));
% I_bic=I+319;


figure;     hold on;    hidden on;
plot(sp(I_fum,:))
hold on;
plot(10*sp(I_mal,:),'color','r')
hold on;
% plot(sp(I_bic,:),'color','g')
% 
% figure;     hold on;    hidden on;
% plot(sp(I_bic,:),'color','b')
% 
 mal2fum=sum(sp(I_mal,:))/sum(sp(I_fum,:))

% bic2pyr=max(sp(I_bic,:))/max(sp(I_pyr,:))
% 
% hoge=sp(I_bic,:);
% 
